{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "62a3bc0d",
   "metadata": {},
   "source": [
    "# Processing of geofeatures for scenarios"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3be71140",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sqlalchemy import create_engine, MetaData\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "from scipy import stats\n",
    "\n",
    "from shapely.geometry import Polygon\n",
    "from shapely.ops import transform\n",
    "import pyproj\n",
    "\n",
    "#from geoalchemy2 import Geometry  # <= not used but must be imported"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "dad031ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3528: SAWarning: Skipped unsupported reflection of expression-based index unique_user_emails\n",
      "  util.warn(\n",
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3528: SAWarning: Skipped unsupported reflection of expression-based index unique_organization_names\n",
      "  util.warn(\n",
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3198: SAWarning: Did not recognize type 'geometry' of column 'extent'\n",
      "  util.warn(\n",
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3536: SAWarning: Predicate of partial index at_most_one_actively_editing_user_per_scenario ignored during reflection\n",
      "  util.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "dict_keys(['public.spatial_ref_sys', 'public.migrations', 'public.users_organizations', 'public.organizations', 'public.users', 'public.roles', 'public.users_projects', 'public.projects', 'public.issued_authn_tokens', 'public.features', 'public.output_results', 'public.scenarios', 'public.users_scenarios', 'public.api_event_kinds', 'public.api_events'])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "api_engine = create_engine(f\"postgres://marxan-api:marxan-api@marxan-postgresql-api:5432/marxan-api\")\n",
    "api_meta = MetaData(schema=\"public\")\n",
    "api_meta.reflect(bind=api_engine)\n",
    "api_meta.tables.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "00c2a9e0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.8/site-packages/sqlalchemy/dialects/postgresql/base.py:3198: SAWarning: Did not recognize type 'geometry' of column 'the_geom'\n",
      "  util.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "dict_keys(['public.spatial_ref_sys', 'public.migrations', 'public.admin_regions', 'public.admin_regions_0', 'public.admin_regions_1', 'public.admin_regions_2', 'public.wdpa', 'public.features_data', 'public.planning_units_geom', 'public.planning_units_geom_square', 'public.planning_units_geom_hexagon', 'public.planning_units_geom_irregular', 'public.scenarios_pu_data', 'public.scenarios_pu_cost_data', 'public.output_results_data', 'public.scenario_features_data'])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "geo_api_engine = create_engine(f\"postgres://marxan-geo-api:marxan-geo-api@marxan-postgresql-geo-api:5432/marxan-geo-api\")\n",
    "geo_api_meta = MetaData(schema=\"public\")\n",
    "geo_api_meta.reflect(bind=geo_api_engine)\n",
    "geo_api_meta.tables.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8c845dd",
   "metadata": {},
   "source": [
    "We have here again 2 steps logic.\n",
    "\n",
    "We need to create the features first, if there is a recipe for such thing.\n",
    "\n",
    "Later we need to attach those features created to our scenario.\n",
    "\n",
    "If a user updates / deletes features from the recepy we need to reflect it and associate them again to the scenario, this steps will need to be triggered.  \n",
    "Also, will we be deleting those from the features table if a user is no longer using it?\n",
    "\n",
    "SPF and FPF are set in the recipe for an scenario.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6885d933",
   "metadata": {},
   "source": [
    "### Recipe definition for feature creation (TBD - with Andrea, Alex and Barre)\n",
    "\n",
    "First approach for recipe definition in [MARXAN 155](https://vizzuality.atlassian.net/browse/MARXAN-155?atlOrigin=eyJpIjoiNjA5YTczNDJkZDU5NGM2OTg4ZDIyZjQ2ZjQ5ZTAxNWYiLCJwIjoiaiJ9)\n",
    "\n",
    "```typescript\n",
    "{\n",
    "  status: 'draft' | 'created'\n",
    "  features: [\n",
    "    {\n",
    "      featureId: string;\n",
    "\n",
    "      /**\n",
    "       * We can either have marxanSettings *or* geoprocessingOperations.\n",
    "       * If marxanSettings is present (*but* we also need to set `kind: 'plain'`\n",
    "       * because of limitations on the way we can describe union types),\n",
    "       * the feature should be used as is (i.e. no splits/intersections),\n",
    "       * with the given settings.\n",
    "       *\n",
    "       * If geoprocessingOperations is present (and `kind: 'withGeoprocessing'`,\n",
    "       * see above), the features generated via geoprocessing operations will be\n",
    "       * used, each with their marxanSettings as specified for each of them.\n",
    "       */\n",
    "      kind: 'plain' | 'withGeoprocessing';\n",
    "      marxanSettings: {\n",
    "        prop: number,\n",
    "        fpf: number,\n",
    "      },\n",
    "\n",
    "      geoprocessingOperations: [\n",
    "        /**\n",
    "         * Either one operation of kind split/v1 or one operation of kind\n",
    "         * stratification/v1.\n",
    "         */\n",
    "        {\n",
    "          kind: 'split/v1',\n",
    "          splitByProperty: string,\n",
    "          splits: [\n",
    "            {\n",
    "              value: string,\n",
    "\n",
    "              marxanSettings: {\n",
    "                prop: number,\n",
    "                fpf: number,\n",
    "              }\n",
    "            }\n",
    "          ]\n",
    "        },\n",
    "        {\n",
    "          kind: 'stratification/v1',\n",
    "          intersectWith: {\n",
    "            featureId: string,\n",
    "          },\n",
    "          splitByProperty: string,\n",
    "          splits: [\n",
    "            {\n",
    "              value: string,\n",
    "\n",
    "              marxanSettings: {\n",
    "                prop: number,\n",
    "                fpf: number,\n",
    "              }\n",
    "            }\n",
    "          ]\n",
    "        },\n",
    "      ]\n",
    "    }\n",
    "  ],\n",
    "};\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a59ef7b",
   "metadata": {},
   "source": [
    "### Feature creation from recipe\n",
    "\n",
    "* We need the recipe and the scenario attached to it.\n",
    "* The recipe status should be set from draft to processing\n",
    "\n",
    "For each feature (as a geofeature in the api feature table) recipe:  \n",
    "1º We need to check if there is already a feature (geofeature in the api) computed with that recepy. If so jump to the next step; attach it to the scenario.  \n",
    "\n",
    "2º If not, creates the new feature and If so jump to the next step; attach it to the scenario.\n",
    "\n",
    "If the status goes from processing to completed what happens  \n",
    "If the status goes from processing to failed what happens  \n",
    "If the status goes from completed to draft what happens  \n",
    "If the status goes from processing to draft what happens  \n",
    "\n",
    "\n",
    "### split rules\n",
    "* Only bioregional tag.\n",
    "* id parent feature.\n",
    "* property name\n",
    "* value to filter by\n",
    "\n",
    "```sql\n",
    "---SPLIT\n",
    "--- insert in features table. link with recipe entities rule feature class name => <operation><parent_name><property><value>, Intersection TBD\n",
    "\n",
    "insert into features ( feature_class_name, property_name, intersection, tag, creation_status, created_by)\n",
    "values ('split_world_eco_dn_181', 'dn', '', 'ecoregion', 'created', 'e9027216-a9f8-4b0f-8fa0-e2d80e2dbab1')\n",
    "\n",
    "\n",
    "--- feature logical reference instead of phisical copy. feature_id --> parent id, we should be able to filter properties like:\n",
    "\n",
    "select * from features_data fd where features_id = 'f01903fe-2e41-423e-8b64-176c77f10be6' and properties @> '{\"dn\":181}';\n",
    "```\n",
    "\n",
    "### Stratification rules\n",
    "* Only bioregional with species tag.\n",
    "* id parent features.\n",
    "* property name\n",
    "* value to filter by\n",
    "\n",
    "```sql\n",
    "--- Stratification: feature id A, feature id B\n",
    "--- insert in features table. link with recipe entities rule feature class name => <operation| stratification><parent_name><property><value>, Intersection TBD\n",
    "insert into features ( feature_class_name, property_name, intersection, tag, creation_status, created_by)\n",
    "values ('stratification_world_eco_dn_181_myspeciesname', 'dn', '', 'ecoregion', 'created', 'e9027216-a9f8-4b0f-8fa0-e2d80e2dbab1');\n",
    "-----\n",
    "insert into features_data (the_geom, properties, source, feature_id)\n",
    "select st_intersection(the_geom, (select the_geom from features_data fd where features_id = '7ce5ee3d-e7ba-4efc-9471-0d9d62a143b4')) as the_geom, \n",
    "(properties||(select properties from features_data fd where features_id = '7ce5ee3d-e7ba-4efc-9471-0d9d62a143b4')) as properties,\n",
    "'intersection',\n",
    "'<new feature id>'\n",
    "from features_data fd \n",
    "where features_id = 'f01903fe-2e41-423e-8b64-176c77f10be6' and properties @> '{\"dn\":181}'\n",
    "and st_intersects(the_geom, (select the_geom from features_data fd where features_id = '7ce5ee3d-e7ba-4efc-9471-0d9d62a143b4'));\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "699ba3b9",
   "metadata": {},
   "source": [
    "### Attach features to our scenario\n",
    "\n",
    "* We need the id of the scenario.\n",
    "* We need the ids of the features to attach to the scenario.\n",
    "* We need the fpf and the target/prop values for each feature\n",
    "\n",
    "Then we need to do:\n",
    "* Link features with the scenario filtering those that are in our study area and set SPF and (target&&prop) (hectares of habitat)\n",
    "* Calculate the area extent of that feature in the area of study.\n",
    "* Intersect them with the PA using filtering from scenario and get the area Protected (should this be against PU protected or through PA?).\n",
    "* Intersect them with the PU? create a view?\n",
    "\n",
    "* With the event logs created by Andrea we should be able to drop the dependency of the created by modify by on the table scenario_features_data\n",
    "\n",
    "@todo: we need to take into account area calculation Andrew uses [EPSG:3410](https://epsg.io/3410) [more info](https://nsidc.org/ease/ease-grid-projection-gt), For tile displaying we will need the data in [EPSG:3857](https://epsg.io/3857), while for data integrity i would keep the projection in [EPSG:4326](https://epsg.io/4326)  I think is a good idea to keep that projection above 4326 or 3857 only at area calculations. Also take into account geometry vs geography types in postgres in order to calculate areas. \n",
    "\n",
    "@todo 2: We need to take into account the difference between (target&&prop) definition (we must have one of them) and the amount. So far on this iteration we are allways understand the target as the portion of the area of the feature to be protected, but from the marxan user manual this can be something completly different. and doesnt need to be consistent between features.\n",
    "```sql\n",
    "-- Scenario data (extension and feature recipe.).\n",
    "select * from scenarios s \n",
    "inner join projects p  on s.project_id = p.id \n",
    "inner join scenarios_feature_recipe f  on f.scenario_id = p.id \n",
    "where s.id = '2a800cc9-b436-4c3d-b781-54b024e3adbb';\n",
    "```\n",
    "```sql\n",
    "-- Create a conexion between Features and scenario.\n",
    "INSERT INTO scenario_features_data (feature_class_id, scenario_id, created_by, total_area) \n",
    "select id, '2a800cc9-b436-4c3d-b781-54b024e3adbb' as scenario_id, 'e9027216-a9f8-4b0f-8fa0-e2d80e2dbab1' as created_by, \n",
    "st_area(st_intersection(ST_GeomFromText('MULTIPOLYGON (((-10 -10, 10 -10, 10 10, -10 -10)))',4326), the_geom)) as total_area, 123 as fpf, 17 as prop  \n",
    "from features_data fd \n",
    "where features_id in ('b9fcd955-7cc9-4c39-9bbe-5047134703d2', 'f01903fe-2e41-423e-8b64-176c77f10be6') \n",
    " and st_intersects(ST_GeomFromText('MULTIPOLYGON (((-10 -10, 10 -10, 10 10, -10 -10)))',4326), the_geom)\n",
    " -- Make this an upsert operation\n",
    " ON CONFLICT (did) DO UPDATE SET dname = EXCLUDED.dname;;\n",
    "```\n",
    "```sql\n",
    "-- Features intersection with wdpa for our scenario.\n",
    "with features_t as (select s.id, p.the_geom from scenario_features_data s\n",
    "            inner join features_data p  on s.feature_class_id = p.id\n",
    "            where s.scenario_id = '84ac0e1e-d5a0-4f58-b3e8-f2409840078f'),\n",
    "features_wdpa as (select features_t.id as feature_scen_id, st_intersection(\n",
    "                st_intersection(\n",
    "                  ST_GeomFromText('MULTIPOLYGON (((-10 -10, 10 -10, 10 10, -10 -10)))',4326),\n",
    "                   features_t.the_geom\n",
    "                  ),\n",
    "                wdpa.the_geom\n",
    "              ) as the_geom, st_area(\n",
    "            st_transform(\n",
    "              st_intersection(\n",
    "                st_intersection(\n",
    "                  ST_GeomFromText('MULTIPOLYGON (((-10 -10, 10 -10, 10 10, -10 -10)))',4326),\n",
    "                   features_t.the_geom\n",
    "                  ),\n",
    "                wdpa.the_geom\n",
    "              ), 3410)\n",
    "            ) as area_protected\n",
    "from features_t\n",
    "left join wdpa\n",
    "  on features_t.the_geom && wdpa.the_geom \n",
    " and st_intersects(ST_GeomFromText('MULTIPOLYGON (((-10 -10, 10 -10, 10 10, -10 -10)))',4326), wdpa.the_geom)\n",
    " )\n",
    " UPDATE scenario_features_data\n",
    "SET (current_pa) = (select current_pa from (select sum(area_protected) current_pa, feature_scen_id from features_wdpa group by feature_scen_id) s where  feature_scen_id = scenario_features_data.id );\n",
    " \n",
    "```\n",
    "```sql\n",
    "-- Features intersection with our scenario planing units. \n",
    "-- This needs to be used for: the tiles and for planningUnitVSConservationFeature.dat file generation taking into account the ammount (area feature in the PU)\n",
    "with f as (select * from scenario_features_data s \n",
    "            inner join features_data p  on s.feature_class_id = p.id \n",
    "            where s.scenario_id = '2a800cc9-b436-4c3d-b781-54b024e3adbb')\n",
    "select *, st_area(st_intersection(st_intersection(ST_GeomFromText('MULTIPOLYGON (((-10 -10, 10 -10, 10 10, -10 -10)))',4326), f.the_geom), spud.the_geom), true) \n",
    "from f\n",
    "left join scenario_pu_data spud on f.the_geom && spud.the_geom;\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63997cc7",
   "metadata": {},
   "source": [
    "Once all of this is done we will be able to generate the next requires files for marxan: \n",
    "* planningUnitVSConservationFeature\n",
    "* conservationFeature\n",
    "* planningUnits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f1e7ddb",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
